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Abstract 

A new and very general technique for simulating solid-fluid suspensions is 
described; its most important feature is that the computational cost scales 
linearly with the number of particles. The method combines Newtonian dy- 
namics of the solid particles with a discretized Boltzmann equation for the 
fluid phase; the many-body hydrodynamic interactions are fully accounted for, 
both in the creeping-flow regime and at higher Reynolds numbers. Brownian 
motion of the solid particles arises spontaneously from stochastic fluctuations 
in the fluid stress tensor, rather than from random forces or displacements 
applied directly to the particles. In this paper, the theoretical foundations 
of the technique are laid out, illustrated by simple analytical and numerical 
examples; in the companion paper, extensive numerical tests of the method, 
for stationary flows, time-dependent flows, and flnite Reynolds number flows, 
are reported. 
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I. INTRODUCTION 



Numerical simulations, which take explicit account of the hydrodynamic forces between 
the suspended particles, are becoming useful tools for studying the dynamical and rheolog- 
ical properties of suspensions. There are at least three important flow regimes which can 
be addressed by numerical simulation: colloidal suspensions of sub-micron sized particles, 
where Brownian forces and viscous forces balance; suspensions of macroscopic particles [i.e. 
larger than fO/xm), where the viscous forces alone are important; and flows at small but non- 
zero Reynolds number (f < i?e < fOO). At present, computational cost is the limiting factor; 
even with supercomputers, it is not feasible to simulate more than about f 00 particles with 
current methods. Thus, development of reliable and more efficient simulation techniques, 
able to cope with thousands of suspended particles, would have a significant impact on our 
understanding of particulate suspensions, complementing present experimental and theoret- 
ical knowledge. In these two papers a new simulation technique for particulate suspensions 
is described; it combines Newtonian dynamics of the solid particles with a discretized Boltz- 
mann model (McNamara and Zanetti, 1988; Higuera et al., 1989) for the fiuid. The basic 
idea is illustrated in Fig. I, which shows 5 solid particles suspended in a background fiuid. 
The fiuid can be modeled as a continuum (Fig. la), as a molecular liquid (Fig. lb), or as a 
discrete velocity (lattice) gas (Fig. Ic). Because of the large scale separations, the dynamics 
of the solid particles are largely independent of the detailed mechanics of the suspending 
fiuid. Discrete velocity models of the fiuid combine most of the features of a fully molecular 
simulation of solid and liquid phases, but are orders of magnitude faster computationally. 
They have many advantages over conventional methods of simulating particulate suspen- 
sions, which are usually based on complicated, computationally intensive solutions of the 
Stokes equations. By contrast, lattice-gas/lattice-Boltzmann simulations are fast, fiexible, 
and simple. The new method is closely related to earlier suspension modeling using lattice- 
gas cellular automata (Ladd et al., 1988; Ladd and Frenkel, 1989; Ladd and Frenkel, 1990; 
Ladd, I99I; van der Hoef et al., I99I), but the large and uncontrollable statistical fiuctu- 
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ations present in lattice-gas models are suppressed, reducing the need for computationally 
expensive ensemble averaging. 

In the classical theory of suspensions (Happel and Brenner, 1986), the hydrodynamic in- 
teractions are assumed to be fully developed; in other words, there is a complete separation 
of time scales between the dynamics of the fluid and the motion of the solid particles. Most 
simulation methods (for instance Ermak and McCammon, 1978; Brady and Bossis, 1988; 
Ladd, 1988; Tran-Cong and Phan-Thien, 1989; Karrila et al., 1989) utilize this approxima- 
tion, even though it imposes a crippling numerical burden, associated with the global nature 
of the interactions; in such cases one must either make further, drastic simplifications, as 
in Ermak and McCammon, 1978, or pay the steep computational cost of an algorithm that 
scales as at least the square and often as the cube of the number of particles (Brady and 
Bossis, 1988; Ladd, 1988; Tran-Cong and Phan-Thien, 1989; Karrila et al., 1989).^'^ In real- 
ity, hydrodynamic interactions develop in time and space from purely /oca/ forces, generated 
at the solid-fiuid surfaces, which then diffuse throughout the fiuid. Recently a new simula- 



^The exact scaling depends on the problem to be solved. In some instances, particle velocities 
are computed for a given set offerees, in which case the computation can scale as iV^. However in 
many cases, for instance to simulate Brownian motion (Bossis and Brady, 1987), the full 6N X 6N 
diffusion coefficient matrix is needed; here the computational cost is of order iV"^. Moreover, 
determining lubrication forces (Durlofsky et al., 1987) also involves an order iV"^ calculation of the 
friction coefficient matrix. 

^The Stokes equations can be solved very efficiently by spectral methods (Fogelson and Peskin, 
1988; Sulsky and Brackbill, 1991). However, there are then the usual difficulties associated with 
incorporating solid boundary conditions into spectral codes; in the referenced work the particle 
surfaces are represented by point forces. Numerical studies of the Navier-Stokes equations indicate 
that, even for pure fluid flows, the lattice-Boltzmann equation is quite competitive with the best 
spectral methods (Chen et al., 1992). 
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tion technique for particulate suspensions has been developed (Ladd, 1993b) which exploits 
this spatial and temporal locality; as a result the computational cost scales linearly with 
the number of particles. The method is also very flexible; the particle size and shape, the 
electrostatic interactions, the flow geometry, the Peclet number (the ratio of viscous forces 
to Brownian forces), and the Reynolds number (the ratio of inertial forces to viscous forces), 
can all be varied independently. 

There is a well established connection between the dynamics of a dilute gas and the 
Navier-Stokes equations (Chapman and Cowling, 1960). Thus, if we can determine the time 
evolution of the one-particle velocity distribution function n(r, v,t), which defines the den- 
sity of particles with velocity v around the space-time point(r,t), then the hydrodynamic 
fields can be readily evaluated. By introducing the assumption of molecular chaos, i.e. that 
successive binary collisions in a dilute gas are uncorrelated, Boltzmann was able to derive 
the integro-differential equation for n named after him (see Chapman and Cowling, 1960). 
The Navier-Stokes equations follow directly from the Boltzmann equation in the limit that 
the dimensions of the macroscopic fiow fields are much larger than the mean free path be- 
tween molecular collisions (Chapman and Cowling, 1960). Because of its complexity, there 
are few direct numerical solutions of the Boltzmann equation (an exception is Yen, 1984), 
but stochastic, particle-based simulations are quite commonly used in molecular gas dynam- 
ics; the merits of this approach are discussed in Bird, 1976 and Bird, 1990. A variant of 
this approach has recently been introduced to simulate particle suspensions (Hoogerbrugge 
and Koelman, 1992; Koelman and Hoogerbrugge, 1993). However, the one-particle distri- 
bution function n(r, v,t) contains much more information than is strictly necessary to solve 
fiuid-dynamics problems; hence, in recent years, there has been a resurgence of interest in 
discrete- velocity or lattice-gas models (Frisch et al., 1986; Frisch et al., 1987), in which the 
continuous distribution of molecular velocities is replaced by a few, carefully chosen, discrete 
values. An important consequence of this work is the realization that a Boltzmann equation, 
fully discretized in coordinate space, velocity space, and time, is equivalent to an explicit, 
second-order, finite-difference approximation to the Navier-Stokes equations (McNamara 
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and Alder, 1992). Numerical studies have shown that a lattice-Boltzmann approximation 
can be comparable in accuracy and computational cost to state-of-the-art Navier-Stokes 
solvers, either finite difference (McNamara and Alder, 1992) or spectral (Chen et al., 1992). 

In the lattice-Boltzmann approximation, the fundamental quantity is the discretized one- 
particle velocity distribution function ni(r,t), which describes the number of particles at a 
particular node of the lattice r, at a time t, with a velocity c^; r, t, and Ci are discrete, 
whereas rii is continuous. The hydrodynamic fields, mass density /?, momentum density 
j = /ju, and momentum fiux 11, are moments of this velocity distribution: 



For simulations of particulate suspensions, the lattice-Boltzmann approximation has two 
important advantages over a finite-difference approximation to the Navier-Stokes equations. 
First, the connection to molecular mechanics makes it possible to derive simple local rules 
for the interactions between the fiuid and the suspended solid particles (Ladd, 1993b); this 
was demonstrated in our earlier lattice-gas simulations (Ladd and Frenkel, 1990; Ladd, 1991; 
van der Hoef et al., 1991). Second, the discrete one-particle distribution function rii contains 
additional information about the dynamics of the fiuid beyond that contained in the Navier- 
Stokes equations; in particular, the fiuid stress tensor, although dynamically coupled to 
the velocity gradient (Frisch et al., 1987), has an independent significance at short times. 
This additional fiexibility allows us to simulate molecular fiuctuations, leading to Brownian 
motion of the suspended particles. To do this, a random fiuctuation, uncorrelated in space 
and time, is added to the fiuid stress tensor (Ladd, 1993b); the variance of the fiuctuations 
defines the effective temperature of the fiuctuating fiuid (Landau and Lifshitz, 1959). This 
approach is quite different from Brownian dynamics (Frmak and McCammon, 1978) or 
Stokesian dynamics (Bossis and Brady, 1987), where random fiuctuations are applied directly 
to the particles; to include hydrodynamic interactions in these methods requires sampling 
from the 6A^ X 6A^ diffusion matrix, which is extremely time consuming. 

The layout of the remainder of this paper is as follows. In section II the lattice-Boltzmann 



nr-, j = XI 
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approximation is described and its connection to Navier-Stokes fluid dynamics is established. 
Section III describes the implementation of the solid-fluid boundary conditions at the mi- 
croscopic level, together with analytic and numerical results for shear flows and channel 
flows. In section IV fluctuations are introduced into the lattice-Boltzmann approximation; 
it is verifled that the simulations satisfy the fluctuation-dissipation theorem and that the 
correct shear viscosity can be obtained from an appropriate Green-Kubo formulae. In the 
companion paper (Ladd, 1993a), referred to hereafter as paper II, results of extensive tests 
of creeping-flow hydrodynamics are reported, for both periodic arrays and random distri- 
butions of spheres; time-dependent and flnite Reynolds number flows are also discussed. 
Three-dimensional simulations of up to 1024 colloidal particles, moving under the action of 
Brownian forces, are also reported. 



II. DISCRETE BOLTZMANN APPROXIMATION 

The computational utility of the lattice-Boltzmann equation is related to the realiza- 
tion that only a small set of discrete velocities are necessary to simulate the Navier-Stokes 
equations (Frisch et al., 1986). It may be helpful in what follows to imagine an underlying 
mechanical model in which identical particles move with discrete velocities from node to 
node of a regular lattice; much of the kinetic theory dilute gases can then be carried over 
directly to the discretized version. The speciflc model used in this work has 18 different 
velocities corresponding to the near-neighbor and second-neighbor directions of a simple 
cubic lattice. Thus there are six velocities of speed 1, corresponding to (100) directions in 
the lattice and 12 velocities of speed \/2, corresponding to the (110) directions, for a total 
of 18. All quantities in this paper are expressed in "lattice units", for which the distance 
between nearest-neighbor nodes and the time for the particles to travel from node to node 
are both unity. Note that the velocities are such that all particles move from node to node 
simultaneously. 

The time evolution of the distribution functions rii is described by a discrete analogue of 
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the Boltzmann equation (Frisch et al., 1987), 

n,(r + c„ t + l)= n,(r, t) + A,(r, t), (2.1) 

where is the change in rii due to instantaneous molecular collisions at the lattice nodes. 
The post collision distribution ni-\-Ai is propagated for one time step, in the direction c^. The 
collision operator A(n) depends on all the n^'s at the node, denoted collectively by n(r;t); 
it can take any form, subject to the constraints of mass and momentum conservation. An 
exact expression for the Boltzmann collision operator has been derived for several different 
lattice-gas models (Frisch et al., 1987; McNamara and Zanetti, 1988), under the usual 
assumption that the distribution functions n(r; t) are uncorrelated from those at previous 
times. However, such collision operators are complex and ill-suited to numerical simulation. 
A computationally useful form for the collision operator can be constructed by linearizing 
the collision operator about the local equilibrium n^' (Higuera et al., 1989), i.e. 

A,(n) = A.^n^") + ^ A-,(n, - ), (2.2) 

i 

where C is the linearized collision operator, and Ai(n^') = by definition. It is not necessary 
to construct a particular collision operator and from this calculate C; rather it is sufficient 
to consider the general principles of conservation and symmetry and then to construct the 
eigenvalues and eigenvectors of C. However, before doing this, we will determine the proper 
form for the equilibrium distribution function. 

To establish the connection between molecular mechanics and fiuid dynamics, it is usual 
to split the distribution function into an equilibrium part and a non-equilibrium part 

= nr + nr. (2.3) 

The equilibrium distribution is a coUisional invariant [i.e. Ai(n^') = 0), and depends only on 
the local hydrodynamic variables (mass density, stream velocity and, in some cases, energy 
density); for a molecular gas n^' is the Maxwell-Boltzmann distribution. It is well known 
that a Maxwell-Boltzmann local equilibrium leads to the Fuler equations of hydrodynam- 
ics (Chapman and Cowling, 1960); however the most-probable (equilibrium) distribution 
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functions of discrete velocity lattice gases give rise to density dependent advection veloci- 
ties and velocity dependent pressures (Friscli et al., 1987). Therefore we seek a constrained 
equilibrium distribution for our discretized velocity model, that will lead to the correct 
macroscopic fluid dynamics at the Euler level; the required form for the moments of the 
velocity distribution function defined in Eq. (1.1) are 

P = Y.<^ j = E<'c., n^^ = ^nrc,c,=pl + ^uu, (2.4) 

i i i 

where p is the pressure and 11^' is the non-dissipative part of the momentum fiux. As in 
the usual kinetic theory of gases, the viscous fiuxes come from the non-equilibrium part of 
the distribution function. 

The equilibrium distribution can be expressed as a series expansion in powers of the fiow 
velocity u, 

nl'^ = p + al'u ■ Ci + al'uu : c^Ci -|- ag u^j , (2-5) 

where uu = uu — (1/3)m^1 is the traceless part of uu; Eq. 2.5 has the same functional form 
as a small u expansion of the Maxwell-Boltzmann distribution function. The moments of 
the distribution function (Eq. 2.4) can be expressed in terms of the coefficients a'^\ which 
are functions only of the speed, q; 

p-^ J2 = 6al + 12a^ + [g^ + 12a/^] u\ (2.6) 

i 

P~' E = [2«1 + 8a/^]u„, (2.7) 

i 

J2 4 = 6al + 24a^ + [g^ + 24a/^] u\ (2.8) 

i 
i 

+ 4a/^ [-^afi-,S - (l/3)(5„^(5^5 + Sa^Sfjs + SaS^fj-/]} U^Ug, (2.9) 

i 

+ iaf^ [-SaPfS + SapS^s + Sa^Sps + SasSp^jjug. (2.10) 

The tensor Sap-ys is unity when all the subscripts are the same [i.e. S^xxx = ^yyyy = S^zzz = 1) 
and zero otherwise; Sap is the Kronecker delta. The third-order moments (Eq. 2.10) do not 
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contribute to the Euler equations, but they must be proportional to the fourth-rank identity 
tensor {SapS^s + ^a-ySps + SaS^p-y)] otherwise the viscous part of the momentum flux will not 
be isotropic (see Eq. 2.28). A comparison of Eqs. 2.6-2.9 with Eq. 2.4, together with the 
isotropy condition Eq. 2.f0, is sufficient to determine all the coefficients, 

al =(2-3c2)/6, al = f/6, a] = f/4, 4 =-f/6, 

= (3c2 - f)/f2, af =1/12, af=l/8, af= 1/12. 

The definition of the pressure (Eq. 2.4) indicates that it is proportional to the density, i.e. 
p = {l/'i){Y,iiT-T'^i ~ P'^'^) — P'^li later it will be shown that is the speed of sound (see 
Eq. 2.26), as in a normal gas. Since the distribution function must always be positive, 
the speed of sound is bounded by the limits f/3 < < 2/3. However, as the sound 
speed approaches either of the two limits, the simulation becomes unstable with respect to 
variations in fiuid velocity. In our simulations, the intermediate value = \Jl/2 is used, to 
maximize the stability with respect to variations in fiow velocity; in this case 

aJ = I/I2, a^=I/24. (2.12) 

Having constructed an equilibrium distribution appropriate for the inviscid (Euler) equa- 
tions, let us consider next how to obtain the correct form for the viscous terms in the fiuid 
equations. We require that the linearized collision operator satisfy the following eigenvalue 
equations; 

= 0, ^ CjXij = 0, ^CjCgXij = AcjCj, ^Cj^£ij = A^c^. (2.13) 

i i i i 

The first two equations follow from conservation of mass and momentum and the last two 
equations describes the isotropic relaxation of the stress tensor; the eigenvalues A and A^ 
are related to the shear and bulk viscosities (Eq. 2.34). Equations 2.13 account for 10 of the 
18 eigenvectors of C. The remaining 8 modes, comprising some higher-order moments of £, 
are not relevant to simulations of the Navier-Stokes equations and will be ignored. Their 
eigenvalues are set to -I so that these modes are then projected out entirely from the post- 
collision distribution; this both simplifies the simulation and ensures the fastest possible 
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relaxation ol the non-hydrodynamic modes. The computational procedure to update the 
lattice-Boltzmann equation is therefore quite straightforward. At each site the moments /?, 
j, and n (Eq. 1.1), and the equilibrium momentum flux 11^' (Eq. 2.4) are calculated. The 
momentum flux is updated according to Eq. 2.13 

= K% + (1 + A)(n„^ - TT:'^) + (1 + \Bm,, - n7^)(i/3)^„^; (2.14) 

then the post-collision distribution, rii + Ai(n), is determined from the inverse of Eq. 1.1, 

n, + Aj(n) = a^p + a^jaC^a + a^Tl'„pC~c~^ + a^{W^^ - 3/jc^). (2.15) 

The term —3pcla'^' keeps 11 orthogonal to p. 

Next we examine the macrodynamical behavior arising from the lattice-Boltzmann equa- 
tion; our method of solution is the usual multi-time-scale analysis (Erisch et al., 1987). We 
begin with conservation equations for the moments of the distribution function 

^n,(r + c„t + l) =^n,(r,t), (2.16) 

i i 

^nj(r -M)q„ =^nj(r,t)Q„, (2.17) 

i i 

^nj(r + c,,t + l)c~c~^= ^nj(r,t)c~c^-F A ^ n"^''(r, t)c~c^. (2.18) 

i i i 

Y ndr + c., t + l)c] n.(r, t)c^ + A^ ^ nr(r, t)cl (2.19) 

i i i 

To find the long-time, long-wavelength dynamics, we introduce a scaling parameter e, 
defined as the ratio of the lattice spacing to a characteristic macroscopic length; the hydro- 
dynamic limit corresponds to e ^ 1. In a molecular gas the appropriate scaling parameter 
is the Knudsen number, the ratio of the mean-free path between collisions to the macro- 
scopic length scale. The parameter e plays a similar role to the Knudsen number in the 
Chapman- Enskog method (Chapman and Cowling, 1960); it is used, first of all, to separate 
the relaxation of the equilibrium and non-equilibrium distributions, (c./. Eq. 2.3) 

= nf + enr. (2.20) 

However, because the lattice spacing and the mean-free path are comparable, there are 
additional contributions to the viscous momentum fiux, which do not appear in the ordinary 
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kinetic theory of gases (see Eq. 2.32). In order to remove discrete lattice artifacts from the 
macroscopic equations, it is convenient to define a macroscopic space scale ri = er, and two 
macroscopic time scales ti = et and ^2 = e^t; this enables a separation to be made between 
the propagation of sound iti) and the diffusion of vorticity (^2) (Frisch et al., 1987). 

Expanding the finite differences, Tiiir + c^, t + 1) — ni(r, t), in Eqs. 2.16-2.19 about r and 
t, and collecting terms that are first order in e we obtain the relaxation on the ti time scale; 

5t.E<' + V«E<'c- = 0' (2-21) 

i i 

E + E = 0, (2.22) 

i i 

^tl XI ^T'^^oC^p + E nfc,aC^pC,^ = ^ XI ^T'^'^t^ + X ^T"^ I '^'^ (2.23) 

i i i i 

the gradient operator refers to derivatives on the macroscopic ri space scale, it i.e. V = 9^. 

The equations for mass and momentum conservation (Eqs. 2.21 and 2.22) can be rewrit- 
ten using Eq. 2.4; 

dt,p + V ■{pn) = 0, (2.24) 

dt, (pu) + V • (pi + puu) = 0, (2.25) 

which are the Euler equations of hydrodynamics. Substituting the equation of state p = pc^ 
and linearizing the Euler equations with respect to 6p and u, we find that, for small density 
fiuctuations, 

dlp = cyp. (2.26) 

Equation 2.26 shows that density fiuctuations relax via the propagation of sound waves, on a 
time scale ti, and therefore decouple from the ^2 time scale evolution of the viscous stresses. 
The time derivative that appears in Eq. 2.23 can be evaluated by using Eqs. 2.24 and 2.25 
to express the time derivatives of p and pu in terms of spatial derivatives; 
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i 

= -V^{pUaUpU^) - cl [UaVpp + UpVaP + V^(/JM^)(5„^] . (2.27) 

The spatial derivative ol the third-order moment can be evaluated directly from Eqs. 2.10 
and 2.11, 

^ n""^ C,aC^pC,^ = (1/3) [VaipUp] + V^(/JM„) + V^(/JM^)(5„^] . (2.28) 

i 

In the incompressible limit, variations in density can be ignored, so that 

9h ^T^zaQp + ^ nf C~C^Q^ = {p/3) [VaUp + V^M„ - (2/3) V^M^] , (2.29) 

i i 

with errors of order Vu"^. Then, from Eq. 2.23, the Navier-Stokes form for the viscous 
stresses can be obtained, 

= - XI K^'^'^u^ = -(/o/3A)(V„M/3 + VpUa - (2/3) V^M^(5„^) ; (2.30) 

i 

hereafter the non-equilibrium (viscous) contributions to the pressure will be ignored. Since 
the Mach number in our simulations is typically 10~^ or less, we can safely ignore the effects 
of compressibility. It is interesting to note that if the speed of sound were set to ^1/3 
instead of ^1/2, then inspection of Eqs. 2.27 and 2.28 indicates that the correct form for the 
viscous stresses, including the non-equilibrium pressure, would be obtained (with corrections 
of order Vu"^), even for non-zero Mach numbers. Eor such simulations, an additional density 
of zero- velocity particles is required to maintain stability (McNamara and Alder, 1993). 

The ^2 relaxation of the mass and momentum densities can be found from the order 
terms in the expansion of Eqs. 2.16 and 2.17; 

5t.E<' = ^t2/> = 0, (2.31) 

i 

9t2 ^fCja + (l/2)V^(9ti Y ^T'^^oC^p + ^ uf C,aC^pC,^) + ^ riY'"^ C,aC^p = 0. (2.32) 

i i i i 

Equation 2.31 shows that the fluid is incompressible on the ^2 time scale; all relaxation of 
density fluctuations takes place on the ti time scale. Using Eqs. 2.29 and 2.30 we can express 
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the long-time variation ol the viscous stresses in a form identical to the incompressible 
Navier- Stokes equations, 

dt,{pu) = riV^u (2.33) 

where 

?7 = -/>(2/A + l)/6 (2.34) 

is the shear viscosity; onece again terms proportional to V • u are neglected. 

The shear viscosity (Eq. 2.34) contains two distinct contributions: the first, proportional 
to A~^, arises from molecular-like collisions (Eq. 2.13); the second term comes from the dif- 
fusion of momentum caused by the finite lattice (Eq. 2.29). In most situations of practical 
interest, the coUisional and lattice contributions to the viscosity are of comparable mag- 
nitude; fortunately they both have exactly the same dependence on velocity gradient, so 
that they may be combined into a single transport coefficient. A linear stability analysis 
shows that A must be bounded in the range —2 < A < (Higuera et al., 1989; McNamara 
and Alder, 1993), otherwise the shear stress grows exponentially in time. The bounds on A 
correspond to the simple physical requirement that the viscosity is positive. Actually the 
viscosity is bounded by a more stringent non-linear stability criterion, which has not yet 
been worked out in detail. Essentially the quantity where u is a characteristic fiow 

velocity, must be larger than some positive constant; this means there must always be some 
small but finite damping of the viscous modes. 

Combining the relaxation of the momentum density on the ti and ^2 time scales, we 
obtain the incompressible Navier-Stokes equation 

dtipu) + W ■ {puu) = -Wp + rjV'^u, (2.35) 

with equation of state p = pc^. Once again we again point out that for the 18 velocity 
model used in this work, the simulations are only valid at low Mach numbers; slightly more 
complex models are needed to capture compressibility effects correctly. In the remainder of 
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this paper, it will be assumed that the simulations will be run under conditions of low Mach 
number, with particle velocities U much less than the sound speed c^; thus V • u = to a 
good approximation. 

Many flows involving particulate suspensions occur at low Reynolds number, and can be 
modeled by the creeping-flow or Stokes equations 

V • u = 0, Wp = TyV^u; (2.36) 

or, in terms of the momentum density j = /ju, 

V-j = 0, Vp = vV^l (2.37) 

In our simulations, we do not model the Stokes equations directly, but rather as a limit (for 
small particle velocities) of the linearized, incompressible Navier-Stokes equation 

5tj = -Vp + //V2j; (2.38) 

Eq. 2.38 can be simulated directly by a change in the equilibrium distribution (c./. Eq. 2.5), 

nf = aSV + «?j-c,. (2.39) 

Einally, a signiflcant simpliflcation of the code occurs when A = — 1, corresponding to a 
viscosity rj = p/6. Although such a large viscosity is not suitable for high Reynolds number 
flows, in the creeping flow limit it allows for a considerable simpliflcation of the collision 
operator 

+ A,(n) = aljV + «?j • c„ (2.40) 

which requires about half the number of floating-point operations as Eq. 2.15; most of our 
i?e = simulations use this viscosity. 

III. SOLID-FLUID BOUNDARY CONDITIONS 

To simulate the hydrodynamic interactions between solid particles in suspension, the 
lattice-Boltzmann approximation must be modifled to incorporate the boundary conditions 
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imposed on the fluid by the sohd particles. The basic methodology is illustrated in Fig. 2. 
The solid particles are deflned by a boundary surface, which can be of any size or shape; 
in Fig. 2 it is a circle. When placed on the lattice, the boundary surface cuts some of 
the links between lattice nodes. The fluid particles moving along these links interact with 
the solid surface at boundary nodes placed halfway along the links. Thus we obtain a 
discrete representation of the particle surface, which becomes more and more precise as the 
particle gets larger. Lattice nodes on either side of the boundary surface are treated in 
an identical fashion, so that fluid flUs the whole volume of space, both inside and outside 
the solid particles. However, because of the relatively small volume inside each particle, 
the interior fluid relaxes quite quickly to rigid-body motion, characterized by the particle 
velocity and angular velocity. Thus, on physically important time scales, the interior fluid 
only contributes an added inertial mass to the solid particle. 

In comparison with our previous work (Ladd et al., 1988; Ladd and Frenkel, 1989; Ladd 
and Frenkel, 1990; Ladd, 1991; van der Hoef et al., 1991), here we have chosen to place 
the boundary nodes on the links connecting the interior and exterior regions, whereas in 
our lattice-gas simulations they were located on the nodes closest to the boundary surface. 
There is little to choose between the two methods; the link method has the advantage that it 
provides a somewhat higher resolution of the solid boundary surface, as can be seen (Fig. 2) 
from the much larger number of boundary nodes compared with the number of lattice nodes 
just inside the surface. On the other hand the node method is faster, although this is of little 
signiflcance in the computationally more intensive lattice-Boltzmann (as opposed to lattice- 
gas) simulations. Although at this point in time, our lattice-Boltzmann simulations have 
been limited to simple symmetrical objects; spheres, disks and plane walls, this restriction 
is not fundamental: in fact a limited number of lattice-gas simulations containing elongated 
objects have already been reported (van der Hoef, 1992). 

At each boundary node there are two incoming distributions ni(r,t_|_) and ^^/(r -|- Ci,t_|_), 
corresponding to velocities Ci and Cj/ (cj/ = —Ci) parallel to the link connecting r and r -|- c^; 
the notation ni(r,t_|_) = ni(r,t) -|- Ai(r,t) is used to indicate the post-collision distribution 

15 



(Eq. 2.15). In some cases boundary nodes for two different link directions, perpendicular to 
one another, may be coincident (see Fig. 2); these are treated independently. The velocity 
of the boundary node U;, is determined by the solid particle velocity U, angular velocity fi, 
and centre of mass R, 

U6 = U + fix(r + |c, -R). (3.1) 

By exchanging population density between rii and n^/ we can modify the local momentum 
density of the fluid to match the velocity of the solid particle surface at the boundary node, 
without affecting either the mass density or the stress, which depend only on the sum ni-\-nii. 
Because the stress tensor is unaffected by the boundary-node interactions, it then follows 
that the hydrodynamic stick boundary condition applies right up to the solid surface, without 
any intervening boundary layer (Ladd and Frenkel, 1990). This point will be illustrated in 
more detail later. The mechanism for the boundary-node interactions is illustrated in Fig. 
3. In Fig. 3a we see the two incoming populations, ni(r,t_|_) and ^^/(r -|- Ci,t_|_), interacting 
with a stationary boundary node. In this case, the populations are simply reflected back in 
the direction they came from (Frisch et al., 1987; Cornubert et al., 1991), so that 

n,{r + c,,t + 1) = n,,{r + c,,t+), and n,,{r,t + 1) = n,{r,t+). (3.2) 

In Figs 3b and 3c the effects of a moving object can be seen. In addition to reflection, 
population density is now transferred across the boundary node, in proportion to the velocity 
of the node u^, 

n,{r + c,,t + 1) = n,f{r + c,,t+) + 2al' pUb ■ 

n,f{r, t + 1) = n,{r, t+) - 2al' pUb ■ c,; (3.3) 

these results are ensemble averages of our earlier boundary-node collision rules for lattice 
gases (Ladd and Frenkel, 1989; Ladd, 1991). Only the velocity component of the boundary 
node along the link direction (cj) is included in the calculation of population transfer; thus 
the outcome in Figs. 3b and 3c is the same. The general form for the boundary node 
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interactions in Eq. 3.3 is determined by the requirement that the local mass density and 
stress tensor are conserved; thus rearrangements of population can only be made among 
pairs of opposite velocities. Furthermore, for stationary nodes we must recover the usual 
"bounce-back" condition (Eq. 3.2). The exact amount of population density transferred [i.e 
the magnitude of the U;, • Ci term) is determined by the requirement that any distribution 
consistent with the boundary-node velocity U;, is stationary with respect to interactions with 
the boundary nodes. It is not obvious that Eq. 3.3 satisfies this condition, but we will verify 
that this is indeed so in the next paragraph. 

Let us now examine boundary-node interactions in more detail. Erom section II we know 
that the distribution function at a node can be written as the sum of equilibrium (Eq. 2.5) 
and non-equilibrium contributions (Eq. 2.30) 

n"^' = —a^'a : c^Ci. (3-4) 

Ignoring higher order gradients [i.e. VVu) and terms proportional to V • u, the coUisional 
stress tensor in Eq. 3.4 can be expressed in terms of velocity gradients (Eq. 2.30), so that 

An"^' = al' pCiCi : Vu. (3-5) 

Then the post-collision distribution ni(r,t_|_) = nj^'(r,t) -|- (1 -|- A)n"^'(r,t) (c./. Eqs. 2.14 
and 2.15) is given, to the same approximation, by 

n,{r,t+) = n,{r,t) + al' pc,c, : Vu(r) = n,f{r,t) + 2al' pu{r) ■ c, + al' pc,c, : Vu(r), (3.6) 

where we have exploited the symmetries in the distribution functions for velocities i and i'. 
If there is a boundary node located at r -|- |ci, then the population ni/(r,t -|- 1) is modified 
according to Eq. 3.3, i.e. 

n,,{r, t + l) = n,,{r, t) + 2aJv [u(r) + (l/2)c, • Vu(r) - U6(r + |c,)] • c,; (3.7) 

thus the distribution is stationary when the fiuid velocity u(r-|- |ci) = u(r) -|- (l/2)ci • Vu(r) 
is equal to the boundary-node velocity u^. 
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To illustrate the action of the boundary nodes more clearly, we consider, as an explicit 
example, planar Couette flow. Figure 4 shows a two-dimensional projection of the lattice- 
Boltzmann model onto the xy -plane; the system is assumed to be time independent, and 
translationally invariant in the y- and 2;-directions. As an idealized model of a solid particle 
surface, two inflnite planes of boundary nodes are set up, at x = and x = L. In the 
fluid between the boundary surfaces (0 < x < Z) there is a uniform velocity gradient 
xUy{x) = 7; outside the boundary planes, the fluid moves with uniform velocity equal to 
the wall velocity. Note that in this example, the lattice nodes are more conveniently set 
at half-integer values of x. The problem then, is to flnd the distribution function for this 
flow geometry that is stationary under the action of the boundary-node microrules, with 
velocities \ih{x = 0) = and \ih{x = L) = ^L. A similar problem, involving mixed stick-slip 
boundary conditions at a stationary wall, has been addressed by Cornubert et al., 1991. 

The expected velocity distribution in a uniform velocity gradient can be constructed from 
the equilibrium distribution for Stokes flow (Eq. 2.39) and the non-equilibrium distribution 
Eq. 3.5, 

^r'' = p^c.^c.y 1 12. (3.8) 

Using the notation of Eig. 4, the velocity distribution function in the fluid (0 < x < Z) can 
be written explicitly as 

no(x) = 4, 

ni{x) = 4, n2{x) = 4, 

mix) = 4:{1 + 2-fx), n4{x) = 4:{1 -2-fx), (3.9) 

nsix) = (1 + 2-fx + 27/A), neix) = (1 - 2-fx + 27/A), 

nr{x) = (1 - 2-fx - 27/A), ns{x) = (1 + 2-fx - 27/A); 

the densities Uq through ^4 have been multiplied by 4 to account for the number of projected 
velocities, and the mass density has, for convenience, been set equal to 24. The velocity 
distribution away from the boundaries can be updated according to the usual time evolution 
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of the lattice-Boltzmann equation. From Eqs. 1.1, 2.14, and 2.15, we can compute the post- 
coUision distribution and then propagate it to the neighboring nodes using Eq. 2.1. The new 
velocity distribution, denoted by n', is 



nQ[x 
n[[x 

n'^[x 



4, 
4, 



n'^{x) 



4(1 + 27x), 

[l + 27(x-l) + 27(l + A)/A], 
[l_27(x-l)-27(l + A)/A], 



= 4(l-27x), (3.10) 
<(x) = [l-27(x + l) + 27(l + A)/A], 
n;,(x) = [l + 27(x + l)-27(l + A)/A]; 

which is identical to the initial distribution (Eq. 3.9), as required. Eor lattice-nodes adjacent 
to the solid-fluid boundaries, the update of some of the populations densities is affected by 
the boundary nodes (Eq. 3.3): explicitly 



n 



' (i) = [1 - 27(1/2)+ 27(1 + A)/A], 



n 



5\2 



[1 + 27(1/2)- 27(1 + A)/A], 



',{L - I) = [1 + 27(Z - i) + 27(1 + A)/A] - 47Z, 
',{L - I) = [1 - 27(Z - i) - 27(1 + A)/A] + 47Z; 



(3.11) 



which, once again, is identical to the initial distribution (Eq. 3.9). Thus the boundary-node 
collision rules generate an exact linear shear flow; this is because they maintain the second 
order accuracy of the pure fluid model. Eurthermore, the velocity distributions outside and 
inside the particle are isolated from one another; thus a sharp change in velocity gradient 
from the inside to the outside the particle surface can be supported, as illustrated in Eig. 4. 

As a result of the boundary-node interactions (Eq. 3.3), forces are exerted on the solid 
particles at the boundary nodes, i.e. 

f(r + + I) = -K-(r + c,,t + 1) - n,f{r,t + 1) - n,{r,t+) + n,f{r + c,,t+)]c, 

= 2[n,{r,t+) - n,,{r + c,,t+) - 2al' pUb ■ c,]c,; (3.12) 

thus momentum is exchanged locally between the fluid and the solid particle, but the com- 
bined momentum of solid and fluid is conserved. The forces and torques on the solid particle 
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are obtained by summing f(r + |ci) and (r + |ci)xf(r + |ci) over all the boundary nodes 
associated with a particular particle. As an example, Eq. 3.12 can be used to compute the 
drag force per unit area on a planar wall adjacent to a steadily shearing fluid. We compute 
the force on one face of each solid boundary surface, assuming that the fluid on the other 
side is moving uniformly with the velocity of the boundary (as shown in Fig. 4) and there- 
fore exerts no force on the wall; this corresponds to replacing one of the distributions in 
Eq. 3.12 by its equilibrium form (Eq. 2.39) with a velocity equal to the wall velocity. Using 
the distributions at just after the molecular collision process (Eq. 3.9 with 1/A replaced 
by (1 + A)/A) we find for the wall forces 

/,(0) = 2[-ne{lt+) + ns{lt+)] = -4(2/A + 1)7 = 777, 

fyiL) = 2[nsiL - I, t+) - nriL - \, t+) - 47Z] = 4(2/A + 1)7 = -r,r, (3.13) 

the last equality follows from summing the coUisional and lattice contributions to the vis- 
cosity (Eq. 2.34), using p = 24. Thus the wall force is computed exactly for linear shear 
fiows. 

As a preliminary application of the method to time-dependent fiows, we consider the 
evolution of the fiow field from an impulsively started fiat plate. The geometry is similar 
to Eig. 4, with the plates being sufficiently far apart that they do not interact over the 
duration of the simulation. We focus on a single plate [x = 0). Initially the system is at 
rest; at t = 0, an impulsive force gives the plate a constant velocity [0, [/, 0]. In this problem 
it is again convenient to define the lattice nodes at half-integer values of x, and for the fiuid 
to reside at the lattice nodes at half-integer values of the time; then the boundary conditions 
at the plate are applied at x = and t = precisely. We compute numerically the evolution 
of the fiow field [0, u{x, 0' 0], created by the diffusion of vorticity into the fiuid, and compare 
with the analytic solutions for the velocity field (Batchelor, 1967) 

u{x,t) = U[l-^{xlVl^t)], (3.14) 

and the force per unit area 
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fit) = -V^MO, t) = -f]U{Trpt)-'^'; (3.15) 

V = J] I p is the kinematic viscosity of the fluid and $ is the error function. The results are 
shown in Fig. 5 (for a viscosity v = 1/6) at several different times. It can be seen that 
there is complete agreement, except at very short times and distances. This simple test 
implies that we can accurately simulate both stationary and time-dependent flows, as will 
be confirmed by further results in paper II. 

Next we consider fiow perpendicular to the wall. The incompressibility condition means 
that there can be no velocity gradients in steady fiow; thus the fiuid velocity and the wall 
velocities are [u, 0, 0]. For this fiow we can further simplify the model in Fig. 4 to only three 
velocity directions; stationary particles, and particles moving in the positive and negative x 
directions. The projected distributions are 

no(x) = 12, 

ni{x) = 6(1 + 2m), n2{x) = 6(l-2u), (3.16) 

Clearly this distribution is stationary with the respect to the evolution of the lattice- 
Boltzmann equation; it is also stationary with respect to the boundary-node update rules 
(Fq. 3.3). However this is not the only stationary distribution that satisfies the boundary 
conditions. A more general form for the distribution function, which exhibits a two time-step 
repeat cycle is 

no{x) = 12, 

ni(x) = 6(1 + 2u + (-l)*+^2u;), n2{x) = 6(1 - 2u - (-l)*+^2u;); (3.17) 

it is straightforward to verify that this distribution also satisfies both the time evolution 
equation and the boundary conditions at the walls. Thus the fiuid has a uniform momentum 
pu and a "staggered" momentum [ — lY'^^pw. Staggered momenta are an artifact of all 
lattice models (Zanetti, 1989); the precise value of the staggered momentum depends on the 
channel width (in this example) and the initial conditions. Although staggered momenta 
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cannot arise spontaneously in the fluid, they can be generated at sohd surfaces, as we see 
in this example. However, it can be shown, by techniques similar to those used above, 
that staggered momentum parallel to the walls is damped by the boundary conditions, so 
that for the plane Couette flow problem, the steady solutions have no staggered momentum 
component. 

In the more complex geometries that are of interest in particulate suspensions, it is 
impossible to analyze the staggered momenta analytically. Numerical results show that 
large oscillations in particle torques can be built up by a feedback mechanism in which 
the staggered momenta are fed by ever increasing angular velocities of the particles. To 
overcome these instabilities we average the force and fluid velocity over two successive time 
steps which effectively cancels out the staggered momentum contribution. In the above 
example, this gives a uniform flow fleld with velocity [u, 0, 0] regardless of the magnitude 
of the staggered momenta. Since the forces at the boundary nodes are generated at the 
half-integer time steps, the smoothly varying force f at the intermediate integer time is 

f (r + ic„ t) = (I/2)[f (r + ic„ t - i) + f (r + ic„ t + i)]. (3.18) 

We can calculate the smooth part of the fluid velocity fleld at half-integer time steps in a 
similar way 

u(r,t + i) = (l/2)[u(r,t) + u(r,t + I)], (3.19) 

or at integer time steps using the 3-point formula 

u(r,t) = (I/4)[u(r,t - I) + 2u(r,t) + u(r,t + I)], (3.20) 

The velocities of flnite mass particles (as opposed to inflnitely massive flxed objects) are 
updated every two time steps, 

U(t + I) = U(t- I) + 2M-^F(t), fi(t + I) = fi(t - I) + 21"^ • T(t). (3.21) 

The particle mass M and moment of inertia I are preassigned parameters which control the 
rate at which particles respond to the fluid flow; usually M and I are on the order of several 
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thousands (in lattice units). Since the velocities vary slowly on the time scale of a lattice- 
Boltzmann cycle, the precise form for the update is usually not too important; however, it 
is important to use time-smoothed forces and torques, as described in Eq. 3.18. 

So far we have considered only linear shear flows. In the next and last example we 
consider stationary two-dimensional channel or Poiseille flow. The geometry is, once again, 
as shown in Fig. 4, with the walls at rest (u;, = 0). The fluid is driven by a pressure gradient, 
which is represented in the simulation by a uniform force density in the fluid. Thus we apply 
a constant increment Ajy to the y momentum at each node so that the pressure gradient 
down the channel is Vj^p = Aj,^. The fluid velocity at each node is measured after half the 
force has been applied; we have found empirically that this prescription gives the fastest 
convergence as a function of system size. The steady flow profiles are compared with the 
analytic solution. 



in Fig. 6. It can be seen that the agreement is very good for channels more than about 9 
lattice spacings wide. Furthermore, we find that the force on the wall is exactly correct at 
steady state, no matter what the channel width. This result follows from the balance between 
forces on the walls and the total force from the pressure gradient. Since the pressure forces 
are distributed equally on the walls it follows that the wall force per unit area / = AjyL/2^ 
which is the correct result for Poiseille fiow. 

In this section we have seen how moving solid boundaries can be incorporated into a 
lattice-Boltzmann simulation, and we have indicated how they function for a few simple 
examples. In paper II we describe numerical results for spherical particles, both periodic 
arrays and random assemblies. Results are compared with known analytic and numerical 
solutions of the creeping-fiow and Navier- Stokes equations. 
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IV. FLUCTUATIONS 



In recent years, it has become increasingly obvious that the the lattice-Boltzmann equa- 
tion is a much better simulation tool for hydrodynamics than lattice gases. However, in 
its normal state the lattice-Boltzmann equation cannot model the molecular fluctuations in 
the solvent that give rise to Brownian motion. Of course in many situations Brownian mo- 
tion is unimportant, but, for suspensions of sub-micron sized particles, it is a fundamental 
component of the dynamics; It has been recently shown (Ladd, 1993b) that fluctuations can 
be incorporated into the lattice-Boltzmann equation, within the framework of fluctuating 
hydrodynamics (Landau and Lifshitz, 1959), by adding a random component to the fluid 
stress tensor. Numerical tests showed that the resulting particle motion, in dilute to con- 
centrated suspensions, closely matches recent experimental results (Zhu et al., 1992; Kao 
et al., 1993), even at very short times where particle inertia plays an important role. In this 
section, we describe the basic theory of fluctuations as it applies to the lattice-Boltzmann 
model; in paper II numerical tests of the method for particulate suspensions of spheres will 
be reported. 

The basic idea of fluctuating hydrodynamics is that, on length scales and time scales 
intermediate between the molecular and the hydrodynamic, thermally-induced fluctuations 
can be reduced to random fluctuations in the fluxes of the conserved variables; i.e. the 
stress tensor and perhaps the heat flux also. Since the fluxes are included explicitly in a 
lattice-Boltzmann simulation, it is plausible that molecular fluctuations can be modeled re- 
alistically on intermediate scales, even though the details of the microscopic interactions are 
different. In the present context, this means that the time evolution of the velocity distri- 
bution includes a stochastic term n'(r,t), representing the thermally-induced fluctuations in 
the stress tensor. 





where n' is chosen so that only its stress moment is non zero (c./. Eq. 3.4), 
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/ c, / 

n-i — a^pCiaCip. 



(4.2) 



These random stress fluctuations are uncorrelated in space and time (Landau and Lifshitz, 
1959) and are sampled Irom a Gaussian distribution, 



the choice ol the variance A serves to deflne the effective temperature ol the fluid. We can 
determine the relationship between A and the temperature via the fluctuation-dissipation 
theorem. 

It is convenient in what follows to work in Fourier space, deflning the Fourier transform 
of the velocity distribution function as 



rev 

where the sum is over all the lattice points in the periodic unit cell. Then the equation for 
momentum conservation (Fq. 2.17) can be written as 



we will attempt to cast this equation into the form of a Langevin equation for the transverse 
(or solenoidal) momentum fluctuations j_L(k, t) = (1 — kk) • j(k, t) (k is the unit vector k/A;). 
Fxpanding the exponential, we have 



jx(k,t + 1) - ji(k,t) + zk • n(k,t + 1) • (1 - kk) - (A;V6)ji(k,t + 1) = 0{k'), (4.6) 



where n(k,t) = ^i(k, t)ciCi is the traceless part of the momentum flux. To obtain a 
Langevin equation, we must divide the momentum flux into a slowly relaxing dissipative 
part and a rapidly varying fluctuating part. In section II, it was shown that on the shorter 
time scale, the stress relaxed to its Navier-Stokes form (Fq. 2.30), whereas the velocity 
gradients varied on a longer time scale. Thus we deflne the fluctuating stress S as 



(4.3) 




(4.4) 




(4.5) 




(4.7) 
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where we have ignored density fluctuations (because they do not couple to the transverse 
momentum flux) and non-linear puu terms (because we are examining small length scales 
where the Reynolds number is negligible). Combining Eqs. 4.6 and 4.7, and summing the 
coUisional (1/3A) and lattice (1/6) contributions to the viscosity (Eq. 2.34), we obtain a 
discrete Langevin equation for the transverse momentum fluctuations, 

jx(k,t + l) -jx(k,t) + //A;2j^(k,t + l) = zk-S(k,t + l)-(l-kk), (4.8) 

with a random force zk • S(k, t + 1) • (1 — kk). The solution of a discrete Langevin equation 
is described in the Appendix; the result in this case is (Eq. A9) 

oo 

2//(jx(k).jx(-k))= Yl (k-S(k,t).(l-kk).S(-k,0)-k). (4.9) 

t= — oo 

The effective temperature of the fluid can be deflned by equating the fluctuations in 
momentum to those in real fluids; for a molecular fluid of N particles of mass m, 

(ji(k) • ji(-k)) = ^ ^ E e-^'^-^-v, • = 2NmkBT = 2pVkBT, (4.10) 

where the factor 2/3 comes from the missing longitudinal fluctuations. In the long- 
wavelength limit we can isotropically average over the directions of k with the result 

lim (k • S(k, t) • (1 - kk) • S(-k, 0) • k) = i (S(t) : S(0)) = 2 (E.,(t)E.,(0)) , (4.11) 

where 

S(t) = lim S(k, t) = Y^ o-(r, t). (4.12) 

rev 

Combining Eqs. 4.9-4.11 we obtain a fluctuation formula for the viscosity 

oo 

riVkeT = (1/2) (E,,(0)E,,(0)) + ^ (E,,(t)E,,(0)) , (4.13) 

t=i 

which is an analogue of the Green-Kubo relation for molecular liquids (Hansen and McDon- 
ald, 1986), discretized in time. The summation in Eq. 4.13 is the equivalent of a Simpson's 
rule approximation to the time integral that appears in the usual Green-Kubo formulae; in 
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paper II this expression is used to calculate suspension viscosities, replacing the fluid stress 
tensor with the combined stress tensor of the solid and fluid phases. 

Let us now consider the time evolution of the stress tensor in the pure fluid. The total 
stress fluctuations S in a volume V are independent of the propagation of population density; 
they vary only because of collisions and random fluctuations. Thus the time evolution of S 
including the random fluctuations can be written as 

S(t + 1) = (1 + A)S(t) + J2 (4-14) 

rev 

Since (<T'(r, t)I](0)) = for all t > 0, the time correlation function of the stress fluctuations 
can be written in terms of the equal time correlation function (c./. Eq. A4), 

(S(t)S(0)) = (l + A)*(S(0)S(0)); (4.15) 

thus from Eqs. 4.13 and 2.34 we can relate the effective temperature to the equal-time stress 
fluctuations, 

pVkBT = 3{j:l). (4.16) 

Lastly, we need to relate the equal-time stress fluctuations to the fluctuations in the random 
stress tensor. This can be done by noting that Eq. 4.14 is a discretized Langevin equation 
for the stress tensor, in which the random forces at different times are uncorrelated (see 
Appendix). Since | (1 -|- A) |< 1, the equal-time correlation function is given by (Eq. A8) 

AV 

i - [i + A) 
Then, from Eqs. 4.16 and 4.17 

A = ipkBT/3)[l - (1 + Xf] = 2r,kBT\\ (4.18) 

which is the fluctuation-dissipation relation for our fluctuating lattice-Boltzmann equation. 
It defines the effective temperature of the fiuid so that the dissipative (Eq. 2.34) and fiuc- 
tuating (Eq. 4.13) expressions for the viscosity are consistent. 
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Some numerical tests of the fluctuating lattice-Boltzmann equation are illustrated in Fig. 
7. Here we compare simulation results for the stress-stress correlation function with theo- 
retical results (Eq. 4.15), for a range of values of A covering kinematic viscosities from about 
10""^ to about 30. For both the over-relaxing collision operators (—2 < A < —1), where the 
stress tensor changes sign at every time step, and for the under-relaxing collision operators 
( — 1 < A < 0) the agreement between theory and simulation is perfect. In general, we will 
choose values of A close to —1, to minimize the relaxation time of the stress fluctuations. 
The special case A = — 1 is particularly useful; here the stress fluctuations decay instanta- 
neously, since only the random part of the stress tensor gets propagated at each time step 
(see Fq. 4.14). In this limit we flnd that the expression for A (Fq. 4.18) reduces to the 
Landau-Lifshitz result (Landau and Lifshitz, 1959) A = 2r]kBT. 



V. CONCLUSIONS 

In this paper the theoretical foundation for a new simulation technique for particu- 
late suspensions has been described. Much of the theory underlying this work is now well 
understood; in particular the macroscopic fluid dynamics arising from the various lattice- 
Boltzmann models, and the modeling of hydrodynamic stick boundary conditions by local 
modiflcations to particle populations. The computational effectiveness of the method will be 
demonstrated in paper II. The inclusion of fluctuations in the stress tensor is a more recent 
development (Ladd, 1993b) and is less well understood. Nevertheless, we have been able 
to derive discrete analogues of the basic equations of fluctuating hydrodynamics; again the 
numerical tests reported in paper II provide convincing numerical evidence of the correctness 
of the approach. 
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FIGURES 



Fig. 1. Microscopic models of a colloidal suspension. The dynamical properties of the 
large particles are insensitive to the detailed motions of the background fluid, so that the 
continuum fluid in Fig. fa can be equally well replaced by either a molecular solvent (Fig. 
fb) or a lattice gas (Fig. fc). Lattice-gas/lattice-Boltzmann simulations are many orders of 
magnitude faster computationally. 

Fig. 2. Location of the boundary nodes for a circular object of radius 2.5 lattice spacings. 
The velocities along links cutting the boundary surface are indicated by arrows. The location 
of the boundary nodes are shown by solid squares and the lattice nodes by solid circles. 

Fig. 3. Population densities before and after a collision with a boundary node. The 
effects of stationary (Fig. 3a) and moving (Figs. 3b and 3c) boundary nodes on the incoming 
populations are shown. The arrows indicate the velocity direction and the lengths of the 
solid lines are proportional to the population densities. The differences between population 
densities are highly exaggerated for clarity. Note that the effects of the moving boundary 
are the same in Figs 3b and 3c, because the velocity component parallel to the link direction 
is the same. 
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Fig. 4. A two-dimensional projection of the lattice-Boltzmann fluid, bounded by plane 
walls. The boundary nodes are shown as squares and the boundary planes [x = and x = L) 
by dashed lines. The circles are the lattice nodes; the set of nodes explicitly considered in 
the text are shown flUed. The fluid between the walls is subjected to a velocity gradient 
by the relative motion of the walls; the fluid outside the walls moves with uniform velocity 
equal to the wall velocity. The labeling of velocity directions used in the text is also shown; 
velocity components in the 2;-direction have been projected onto the xy-plane. There is an 
additional density of stationary particles (not shown), labeled 0, corresponding to velocities 
[0,0, ±1]. The lower portion of the flgure is a plot of the velocity proflle, with the crosses 
showing the fluid velocity Uy at the nodes and the dotted line the interpolation between 
nodes. 

Fig. 5. Flow induced by an impulsively loaded flat plate, U{t) = for t < and 
U{t) = U for t > 0. The plots show the velocity held and force per unit area on the wall 
at various times; the solid circles are the numerical simulations, and the solid lines are the 
analytical results (Batchelor, 1967). 

Fig. 6. Poiseille flow in two-dimensional channels for different channel widths, L. 
The simulated velocities (solid circles) are scaled by the theoretical maximum value Uq = 

Fig. 7. Stress-stress correlation function for a lattice-Boltzmann fluid. Numerical results 
(solid circles) and theoretical results from Fq. 4.15 (solid lines) are shown for different values 
of A. 
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APPENDIX: DISCRETE LANGEVIN EQUATION 

The discrete Langevin equation (Eqs. 4.8 or 4.14) is of the general form 

jit + l)-jit) = -ajit) + fit), (Al) 

where a is a positive constant controlling the rate of dissipation (0 < a < 2), and f[t) is 
the random force. The random force has the usual property 

imjit')) = 0, (A2) 
for all t > t' . We can rewrite Eq. Al as 

J{t + f) = (1 - ay'j{t) + i:(l - ay-'f{t + t'- s), (A3) 

and from Eq. A2 we can express the time correlation function in terms of equal-time fluc- 
tuations; 

{jit + t')jit)) = il-af{jit)jit)). (A4) 
The equal time correlations, measured from some initial time (t = 0), are given by 

{jit')jit')) = (1 - af o-(o)j (0)) + i: i: (1 - ar-'-' if if - s)fit' - . (A5) 

s=l s' = l 



In the long-time [t' — )• oo) limit the system loses all memory of its initial conditions; the 
equal-time fluctuations can then be written as 

oo oo 

^•') = EE(i-«r^^'(/(^)/(^')); (A6) 

s=0 s'=0 

Equation A6 can be simplified by a change of variables, sj- = s it s'; 

oo oo 

= E (i-«)i^-i(/(3-)/(o)) 

s_ = — OO s-|-=0 

1 oo 

= i E (l-«)l^-l(/(3-)/(0)). (A7) 

Za — 

s_ = — oo 

In the special case that the random force is delta-function correlated in time (as in 
Eq. 4.14), then 
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In the more general case [i.e. in Eq. 4.8), the random force has a finite, although short, 
relaxation time. Here we can only obtain a simple result for small values of a, corresponding 
to the k ^ limit in Eq. 4.8; 

1 CO 

(f) = ^ E imfm , (A9) 

t= — oo 

which agrees with Eq. A8 when a is small and the random force is delta-function correlated. 
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